{
 "metadata": {
  "name": "symmetry"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "raw",
     "metadata": {},
     "source": [
      "This file is an ipython notebook.  Use a browser to edit the file or run any of the code by issuing\n",
      "\n",
      "->ipython notebook symmetry.ipynb\n",
      "\n",
      "from a command line."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# This block builds the linear operator A.  It uses python code first.py and cython code from first_c.pyc\n",
      "# With the parameters specified, the block takes 15 seconds to run.\n",
      "import numpy as np\n",
      "from time import time\n",
      "from first_c import LO\n",
      "u = 2.0e-5  # log fractional error\n",
      "dy = 8.0e-5 # Step size in log of volume     (1.0e-4)\n",
      "n_g = 400   # Number of discrete steps in g (200)\n",
      "n_h = 400   # Number of steps in h           (300)\n",
      "t0 = time()\n",
      "A = LO(u, dy, n_g, n_h) # Integral operator implemented as scipy.sparse.linalg.LinearOperator\n",
      "t0 = time() - t0\n",
      "print('%6.2f seconds to build LO'%(t0,))\n",
      "class LO_B(LO):\n",
      "    def __init__(self, A):\n",
      "        self.A = A\n",
      "        self.shape = A.shape\n",
      "    def matvec(self, v):\n",
      "        A = self.A\n",
      "        return A.rmatvec(A.matvec(v))\n",
      "    def rmatvec(self, v):\n",
      "        A = self.A\n",
      "        return A.matvec(A.rmatvec(v))\n",
      "ATA = LO_B(A)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "n_states=106316 n_pairs=212082029\n",
        " 21.45 seconds to build LO\n"
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Create v, a vector of states with a few selected components set to one and all others set to zero.\n",
      "# Then calculate A(v), A^T(A(v)), and A^T(v).\n",
      "# Fixme: h is mapped to -h in images\n",
      "v = np.zeros(A.n_states)\n",
      "Big = 0.95\n",
      "f_sources = ((-0.98, Big), (-0.65, Big), (-.35, -Big), (0,0), (.5, Big*Big), (0.995, 0.0))\n",
      "i_sources = []                         # For tagging sources in plot\n",
      "for g_,h_ in f_sources:\n",
      "    g = u * g_\n",
      "    h_max = np.sqrt(24*(u-g))\n",
      "    h = h_max * h_\n",
      "    G = A.g2G(g)[0]\n",
      "    H = A.h2H(h)[0]\n",
      "    #print('g_=%5.2f, g=%10.3e G=%3d, h_=%5.2f, h=%10.e, H=%3d'%(g_, g, G, h_, h, H))\n",
      "    i_sources.append((G, H))\n",
      "    k = A.state_dict[(G,H)][-1]        # get 1-d index of state vector that corresponds to (g, h)\n",
      "    v[k] = 1                           # Set component to 1.0\n",
      "t1 = time()\n",
      "f_v = A.matvec(v)                      # Map v forward\n",
      "fb_v = A.rmatvec(f_v)                  # Map f_v backward\n",
      "b_v = A.rmatvec(v)                     # Map v backward\n",
      "t2 = time()\n",
      "print('%6.2f seconds to map forward and backwards.'%(t2-t1,))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "  4.16 seconds to map forward and backwards.\n"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot v, A(v), A^T(A(v)), and A^T(v).\n",
      "import matplotlib as mpl\n",
      "import matplotlib.pyplot as plt\n",
      "z = A.vec2z(np.ones((A.n_states,))) # Make 2-d array of ones and zeros for masking plots\n",
      "def two_d(w):\n",
      "    'return 2-d version of state vector suitable for plotting'\n",
      "    u = A.vec2z(w.reshape((A.n_states,)))\n",
      "    m = u.max()\n",
      "    w = u*z + m*z\n",
      "    for G,H in i_sources:\n",
      "        t = w[G,H]\n",
      "        w[G-1:G+2,H-2:H+3] = 0 # Make big markers for source points\n",
      "        w[G,H] = t\n",
      "    return w\n",
      "fig = plt.figure(figsize=(24,12))\n",
      "fig.suptitle('Integral operator A for dy=%f'%dy)\n",
      "h_max = (48*u)**.5\n",
      "def plot(fig, location, data, title=''):\n",
      "    data = two_d(data)\n",
      "    ax = fig.add_subplot(*location)\n",
      "    ax.imshow(data.T[-1::-1,:], interpolation=\"nearest\", extent=[-u,u,-h_max,h_max], aspect='auto')\n",
      "    ax.set_title(title)\n",
      "    ax.set_xlabel('$g$')\n",
      "    ax.set_ylabel('$h$')\n",
      "plot(fig, (2,2,1), np.ones((A.n_states,1)), '$v$')\n",
      "plot(fig, (2,2,2), f_v, '$A(v)$')\n",
      "#plot(fig, (2,2,3), fb_i, '$A^T(A(v))$')\n",
      "plot(fig, (2,2,4), b_v, '$A^T(v)$')\n",
      "plt.show()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Calculate eigenfunction\n",
      "tol = 1e-6\n",
      "maxiter = 10000\n",
      "t1 = time()\n",
      "w,b = A.power(n_iter=maxiter, small=tol)\n",
      "t2 = time()\n",
      "print('Calculate eigenfunction of A in %6.2f seconds.'%(t2-t1,))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Calculate eigenfunction of A in 146.35 seconds.\n"
       ]
      }
     ],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# Plot eigenfunctions\n",
      "m_g = 50 # Number of g values to plot\n",
      "m_h = 50 # Number of h values to plot\n",
      "f_factor = 1e-45\n",
      "from mpl_toolkits.mplot3d import Axes3D  # Mysteriously necessary for \"projection='3d'\".\n",
      "\n",
      "# Make mesh grid for plots\n",
      "g,h = A.gh(m_g, m_h)\n",
      "G, H = np.meshgrid(g, h)\n",
      "\n",
      "def b2z(b):\n",
      "    '''Calculate z values to make 2-d plot of vector b\n",
      "    '''\n",
      "    if b.min() + b.max() < 0:\n",
      "        b *= -1\n",
      "    b = b/np.sqrt(A.g_step*A.h_step)\n",
      "    floor = f_factor*b.max()\n",
      "    b = np.log10(np.fmax(b,floor))\n",
      "    b = b.reshape((-1,1))\n",
      "    return A.vec2z(b, g, h, floor=np.log10(floor))\n",
      "\n",
      "def plot(z, w, ax, text=''):\n",
      "    surf = ax.plot_surface(G, H, z.T, rstride=1, cstride=1, cmap=mpl.cm.jet, linewidth=1, antialiased=False)\n",
      "    ax.set_xlabel(r'$g$')\n",
      "    ax.set_ylabel(r'$h$')\n",
      "    ax.set_title(r'%s$n_g=%d$ $n_h=%d$ $dy=%8.3g$ $\\lambda=%8.3g$'%(text,n_g, n_h, dy, np.real(w)))\n",
      "    \n",
      "\n",
      "fig = plt.figure(figsize=(24,12))\n",
      "ax1 = fig.add_subplot(1,2,1,projection='3d', elev=15,azim=-45)\n",
      "z = b2z(b)\n",
      "plot(z, w, ax1)\n",
      "ax2 = fig.add_subplot(1,2,2,projection='3d', elev=15,azim=-135)\n",
      "plot(z, w, ax2)\n",
      "plt.show()\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 5
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}